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We introduce a simple stochastic system able to generate anomalous diffusion both for position and 
velocity. The model represents a viable description of the Fermi's acceleration mechanism and it is 
amenable to analytical treatment through a linear Boltzmann equation. The asymptotic probability 
distribution functions (PDF) for velocity and position are explicitly derived. The diffusion process 
is highly non-Gaussian and the time growth of moments is characterized by only two exponents 
and Vv. The diffusion process is anomalous (non Gaussian) but with a defined scaling properties i.e. 
P(|r|,t) = l/t^^F^^drl/r^) and similarly for velocity. 
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About half a century ago, Fermi proposed an accelera- 
tion mechanism for interstellar particles, now referred as 
Fermi's acceleration to explain the very high energy 
of cosmic rays. Nowadays, Fermi's acceleration remains 
one of the relevant explanations for several phenomena 
in plasma physics and astrophysics In Fermi's 

mechanism, at variance with diffusion in real space {e.g. 
Lorentz's gas 0,0)7 it is the velocity that undergoes dif- 
fusion due to the presence of a stochastic acceleration. 
This original idea Jj found also a number of applica- 
tions in the theory of dynamical systems, because it of- 
fers a sirople but not trivial way to generate chaotic sys- 
tems H %Ml 113 mill 113 to describe the dynamics 
of comets jl4| and the motion of charged particles into 
electromagnetic fields [T5l| . 

This letter studies the Fermi's mechanism through a 
simple model, in which a particle may absorb kinetic en- 
ergy (accelerate) through collisions against moving scat- 
terers. The purpose of the model is to provide a de- 
scription of acceleration mechanism without specific as- 
sumptions on the interstellar media (e.g. presence of 
turbulent magnetic fiels or matter fractal distribution). 
Since we are interested mainly in the general features 
of the diffusion process in phase space, we consider only 
non-relativistic classical particles, neglecting the details 
of their interactions with magnetic fields. Our model 
is able to produce an acceleration process characterized 
by anomalous diffusion both in velocity and position: 
(|v(t)|2) ^ ^2.., (|x(i)P) ^ ^2.^ {v, > 1/2, J., > 1/2) and 
non Gaussian distributions. For simplicity we restrict to 
a two dimensional system (x(t) e R^) and circular ob- 
stacles with radius a. These obstacles represent the re- 
gions where the magnetic field irregularities, responsible 
for the scattering, are localized. We shall assume that 
they are randomly distributed in the plane with uniform 
density p. Their velocity V = {V cos(j),V sincj)) is cho- 
sen according to an isotropic probability G{V)dVd(j)/2Tr, 



with the only requirement that {V'^) < oo. We further- 
more assume, that the influence of the particles on the 
magnetic scatterers is so negligible, that the velocity of 
the latter is kept unchanged during collisions. In the ref- 
erence frame of the scatterers, where the energy of the 
particle remains constant, the collision is completely de- 
termined by the trajectory deflection angle and collision 
time. Such a deflection angle a = a(V, v, b, H) depends 
on particle and obstacle velocity, on the impact parame- 
ter b and on the shape of the magnetic flcld H inside the 
obstacle. 

After a collision, the new particle velocity is = 
M(a)vg , where the superscripts A and B stand for "af- 
ter" and "before" the collision respectively, while sub- 
script S refers to velocities in the scatterers reference 
frame and the rotation matrix ]V[(a) describes the scat- 
tering deflection. In the flxed reference frame, v = 
V5 -I- V, we have: 

= V + M(a)(v'^- V). 

A further simplification occurs when the density of the 
scatterers p is small. Indeed, in the low density limit 
{p ^ Iq^ <C a~2), the mean free path is much larger 
than the average distance between nearest obstacles Iq 
and two successive collisions can be safely assumed to be 
independent. Such a simplification allows describing the 
evolution of the system as a stochastic map: 



v„+i = V + M(a)(v„-V) 



tn + Atr, 



(1) 



where At„ — £n/\'Vn\ and in are the time and the 
free path between two consecutive collisions, respectively. 
Under the hypothesis of independent collisions, well veri- 
fied at low densities, a and Ai„ become independent ran- 
dom variables, whose distribution will be determined by 
the distribution of the obstacle velocities V and impact 



parameters b. For sake of simplicity, we shall suppose 
that the probability distribution Q{a) for a is indepen- 
dent of V and is an even function (for symmetry reasons). 

Let us now specify, the probability distribution of At. 
If the scatterers were at rest, a particle with speed v 
would encounter uniformly distributed obstacles, then 
the variable At would follow an exponential distribution 
law: P{At) = Aexp(— AAi) with a decay rate A = 2apv. 
However, in the case of moving obstacles, a particle will 
encounter preferentially those obstacles with a velocity 
direction opposite to its own velocity. Then, the rate A 
has to be replaced by A(v, V) = 2ap\v — V|. Therefore 
the probability that a particle with velocity v undergoes 
a collision with an obstacle of velocity V, after a time At 
from the previous collision is: 

Pv(Ai|v) = 2ap|v - V\G{V) exp{-A(v) At} , (2) 

with, A(v) — 2ap(|v — V|)v where (....)v indicates the 
average on V 16]. Of course more general laws can be 
considered to take into account the presence of a minimal 
collision time, but our simulations show that the results 
are fairly independent of the details of P{At). 

The above system might look somehow trivial when 
considering its statistical properties as function of the 
number of collisions. Noting that v„_|_i is obtained from 
v„ by a random rotation plus a random shift, one can 
expect, from the central limit theorem, that (|v„p) ^ n 
and analogously (|r„p) ~ n, of course the PDF of v„ 
and „ are Gaussian. However, the proper question is 
about the dependence in time rather than in the n, so the 
particles at the same time t experience a different number 
of collisions and therefore the PDF of v(i) and r{t) are 
not expected to be Gaussian. The time after n collisions, 
t = is the sum of random quantities and its 

dependence on n can be worked out via the following 
scaling argument 0|. Since Atk ~ £/ a/ (|vfcp) ^ l/Vk, 
then t ^ J2k=i k~^^'^ ~ y/n; (|r„p) ~ n correspond to 
the scaling laws: 

(Wm-t' and (|r(t)p)^<^ (3) 

In the following, we will show that this "mean field" scal- 
ing behaviour is in fact correct. 

Let us now present a simple but relevant remark: an 
easy computation gives 

(|x(i)-x(0)|2) ~2 f f \\v,{t')\^)R,[t',T)dt'dT (4) 
Jo Jo 

•w\thRx{t,T) — {vx{t)vx{t + T)l {\vx{t)\'^). In Lorentz-gas 
with fixed obstacles and elastic collisions, {\vx{t')\'^) re- 
mains bounded, thus the integration over t' yields a pro- 
portionality to t. In addition stationarity makes Ri{t' , t) 
independent of t' and an anomalous diffusion can orig- 
inate only by long time tail of -R,;(t) (i.e. slower than 
0(t-1)). 

By contrast, in systems with Fermi's acceleration 
which are non stationary, Rx{t, t) depends also on t, and 



the scaling laws (|v(t)|2) ~ t"^"- and (|r(i)|2) ~ t'^"'' , 
could not be trivially related (see 01 . In this case, one 
can only derive the bound: 

Vx <i^v + ^ , 

the equality holding in presence of a very strong time 
correlation among the velocities. 

A Boltzmann-like description of our model can be car- 
ried out under the hypothesis of the independence of con- 
secutive collisions, already employed to write the dynam- 
ics (|1I2|I . The spatial homogeneous probability f{v,t), 
that a particle has a velocity v at time t, evolves accord- 
ing to the equation 

5J(v,i) =-A(v)/(v,t)+ J duT(v,u)/(u,i), (5) 

where, from Eqs. and (O, the transition probability 
T reads 

r(v, u) = (A(u, V)5[v - V - M(a)(u - V)]> v,a • (6) 

Equation at variance with the usual Boltzmann equa- 
tion, is linear because of the independence between ob- 
stacles and particles. Moreover, the isotropy of the 
scatterer distribution insures that T is only a function 
of the moduli v,v^ and of the angle 9 — 9^ between 
V and v^: T = T{v,v^,9 - 9^). Therefore it is 
convenient to expand / in Fourier series of the angu- 
lar variable: f{v,9) = '}2,kfk{v)eyi\){ik9). The hnear- 
ity of eq. © decouples Fourier modes fk and we oh- 
Uin dtfk{v,t) = -X{v)fk+ J dv^Tk{v,v^)fk{v'',t) with 
Tk{v,v^) = Jd9 exp{ik9)T{v,v^,9). 

We can restrict to the asymptotic time behavior be- 
cause the particles accelerate and their distribution is 
rapidly dominated by velocities v ^ V. Then an asymp- 
totic expansion of the operator T is possible in the small 
parameter V/v. 

We first consider the evolution of an isotropic density 
f{v, 9) = f{v) by introducing the PDF of the velocity 
modulus, V — |v|, g{v,t) = 27rz;/(|v|). We can then 
substitute J[v - V - M(a)(v^ - V)] with S[v^ - v - 
(1 — Af(a))V] in Eq. © and, at second order in v/V, we 
obtain ^18] : 

dtgiv,t) = Dvd^.g , (7) 

where D = apaF{V'^)yr and ctf = (1 — cosa)^. The 
asymptotic solution of this Fokker-Planck equation may 
be explicitly derived [Tsj and it converges toward the scal- 
ing function ga{v,t) = h[v / (Dt)] / (Dt) . A direct sub- 
stitution into Eq. ^ shows that the scaling function 
h verifies the equation: + + h{£,) = 0, 

whose solution is h{^) = ^cxp(— ^). Thus ga{v,t) = 
v/{Dt)^ cxp{—v/{Dt)} and the moments of v are given 
by (w") - (n+ !)!(£>*)". 

We now consider the evolution of a non isotropic 
Fourier mode f{v,9) — exp{ik9)fk{v). We can 
then substitute S[v - M(a)v^ - (1 - M(q!))V] with 
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exp{-ika)S[v^ - v - (1 - M{a))'V], in Eq. ® and per- 
forming the expansion at small v/V we obtain 



dtfkiv) = -Lkfk{v) 



(8) 



at leading order in V/v, with = 2apaF,k and (Jp,k — 
(1 — cos(fca))Q (we have exploited even symmetry of the 
a distribution Q). Thus for any initial distribution func- 
tion, fk relaxes exponentially. This is physically clear, as 
few collisions randomize the phase of the velocity. The 
proportionality to v of the relaxation rate is a conse- 
quence of the fact that the number of collisions per unit 
time is proportional to v. We note that the neglected 
higher order terms in V/v, in Eq. (|SJ|, are a transport 
and a diffusion term, not relevant for large v. 

This result shows that any distribution rapidly be- 
comes isotropic, it also allows to compute the velocity 
autocorrelation function and some properties of the spa- 
tial diffusion. For large t, the velocity autocorrelation can 
be expressed in terms of the asymptotic distribution ga 
and of the solution of the Boltzmann equation Q , with 
a delta- like initial condition 5{-v' — v). Using Eq. ©, we 
can compute this solution. Using this result, we obtain 
the expression for the autocorrelation function: 



{n.{t)v^{t + r)) = i{Dtf/{l + DLitrY 



(9) 



The diffusion properties of velocity can be derived di- 
rectly by Eq. ^ computed at i = and from the sim- 
metry x <-> y 



{Mt)\^)^Q{Dtf 



(10) 



The algebraic decay of the function ij^ is fast enough to 
make the integral Q convergent then 



{\x{t)-x{Qt 



2D 
i 

Li 



(11) 



The above results (|9I1UI11|I agree with the simple argu- 
ment leading to Eq. ||2Jl. In order to study the statistical 
properties of the particle position, we consider the evolu- 
tion of the PDF, /(r,v, i), for the velocity and position. 
The generalization of Eq. (O is the inhomogeneous Boltz- 
mann equation where dtf is replaced by dtf{v,v^t) + 
v.Vr/. We introduce the distribution g{r,(j),v,9,t) = 
47r^rw/(r, V, i) in polar coordinates, r = (r cos c^, r cos 0) 
V — {v cos 9, V cos 9). We limit the discussion to the 
isotropic case, that is g{r, 4>, v, 9) — g{r, v,9 — 4>). Again, 
it is convenient to develop the angular dependence in 
Fourier modes: g[r, cf), v, 9) — J2k 9k{r, v) exTp{ik{9 — </))). 

The evolution equation for each mode gk can be writ- 
ten down and studied in the long time limit, but we do 
not report here the detailed derivation. We can prove 
that go is dominant for large t and converges toward a 
scaling function ho[Ly^r/{D^/H),v/iDt)]Ll^^ /{D'^/H^) 
[l^ which verifies the equation 



d^ho 




(12) 



FIG. 1: PDF of rescaled particle velocity u = |v|/^ {w(i)^) 
at times t = 80,320,1280,5120, corresponding to circles, 
squares, diamonds and triangles respectively. The perfect 
collapse is in accordance with the theoretical results. Inset: 
{v{t)"y^" vs t, for n = 2,4,6 from bottom to top; straight 
lines have slope 1 



We were not able to explicitly solve Eq. 112|l . how- 
ever we computed all the moments of hg by recurrence. 
Indeed, if we define ak^n = /q drdvv^r'^ho{r^v)^ we ob- 
tain: {n + fc)afe^„ = k{k + l)afe„i^„ -|- / 4:ak+i,n-2- Us- 
ing that driho{ri,S,) — h{^) = ^exp(— ^), we have 
Ofc.o = (fc -|- 1)!. The recurrence relation then, for in- 
stance, gives: (r^) = 2{\x{t) - a;(0)p) - Dt'^/Li, 
(r*) ~ 8DH^/{3Ll) and so on. 

Numerical simulations confirm our theoretical results 
obtained via the Boltzmann equation for the scaling of 
the moments of position and velocity. Simulations were 
run for a set of iV = 10^ particles starting from ran- 
dom initial positions and velocities well localized in a 
region of the phase space with size small compared with 
Iq and Vq. For a particle with velocity v, we randomly 
select a scatter velocity V from G{V) ~ S{V — Vo)/2n 
and a scattering angle a uniformely distributed in [0, 27r]. 
Then the collision time Atn was drawn from the law Q 
and finally particle position and velocity were updated 
via the rule Q|. Figs. Q and El show, for different 
times, the rescaled PDF for the velocity modulus v = |v| 
and x-coordinate. The collapse of the curves is very 
good and the exponential tails agree with our theoreti- 
cal predictions. The qualitative scenario corresponding 
io = — 1 and exponential tails remains unaffected 
by using generic form of Pv(At, v) provided it keeps the 
exponential decay for large At. 

We stress that the diffusion process is anomalous, since 
Vx is different from 1/2 (we discuss only the diffusion for 
the position but similar considerations hold for veloc- 
ity). The existence of a unique exponent characterizing 
the growth of the different moments (|r(t)|") ~ is 
somehow peculiar, and in ref. |20|| is called weak anoma- 
lous diffusion, because the most general anomalous diffu- 
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FIG. 2: PDF of the rescaled coordinate z = a:/^/{5(tp) at 
times t = 80, 320, 1280, 5120 with the same symbols of fig.[Tl 
Inset: (a:(t)">^/" vs t, for n = 2,4,6 from bottom to top; 
straight lines have slope 1. 

sion (strong anomalous diffusion) implies that (|r(t)|") ^ 
fnuj;{n) T^[^]j Vx{n) a non constant parameter, therefore all 
the PDF's cannot be rescaled onto a single curve. 



The model can be made more realistic taking into ac- 
count possible energy losses due to interactions with the 
medium and energy irradiation. This dissipation is mim- 
icked by rescaling the velocity components by a factor 
^1 — 7 after each collision, where 7 — — \-^\/[tc{v^ + 
v^)] with Tc a typical time of the collision. The presence 
of dissipation introduces neither substantial modification 
on tailed structure of position and velocity PDFs, nor 
change their scaling behavior. Even different forms of 7 
do not affect the global scenario. 

In summary, we introduced a simplified treatable 
model of Fermi's acceleration. It is remarkable that the 
system presents an anomalous (super) diffusion both for 
position and velocity, which is robust under changes of 
the details and it is characterized by only two expo- 
nents Vx — — ^ and by the PDFs' scaling behavior 
P{\v\,t) = l/t-^Fx{\v\/t-''), P(|v|,t) = l/r-'F.dvIA'^"). 
Furthermore, as a consequence of non trivial correlations, 
the two exponents and are not related by a simple 
dimensional argument. 

A.V. and F.C. acknowledge the financial support by 
Cofin Miur 2001 Fisica Statistica di Sistemi Classici e 
Quantistici. F.B. is supported by E.U. Network Stirring 
and Mixing (RTN2-2001-00285). 



[1] E. Fermi, Phys. Rev. 75, 1169 (1949). 

[2] G. Michalek, M. Ostrowski and R. Schlickeiser, Solar 
Physics 184, 339 (1999). 

[3] M. Ostrowski and G. Siemieniec-Ozieblo Astroparticke 
Physics 6, 271 (1997). 

[4] T.J. Freeman and G.K. Parks Journal of Geophysical Re- 
search 105, 15715 (2000). 

[5] H.A. Lorentz, Proc. R. Acad. Sci. Amsterdam 7, 438 
(1905). 

[6] J.R. Dorfman An Introduction to Chaos in Non- 
equilibrium Statistical Mechanics (Cambridge University 
Press; Cambridge UK 1999). 

[7] G.M. Zaslavsky and B.V. Chirikov, Sov. Phys. Dokl. 9, 
989 (1965). 

[8] A.J. Lichtenberg and M.A. Lieberman, Regular and 
Chaotic Dynamics (2-nd Ed. (Springer- Verlag, New 
York, 1991) 

[9] T. Geisel, J. Wagenhuber, P. Niebauer, and G. Obermair, 

Phys. Rev. Lett. 64, 1581 (1990). 
[10] V.V. Afanasiev, R.Z. Sagdeev, and G.M. Zaslavsky, 

Chaos 1, 143 (1991). 
[11] N.D. Armstead, B.R. Hunt and E. Ott Phys. Rev. E 67, 

021110 (2003). 

[12] M.Y. Natenzon, A.I. Neishtadt, R.Z. Sagdeev, 
G.K. Seryakov, G.M. Zaslavsky, Phys. Lett. A 145, 255 
(1990). 

[13] H. Kunz, R. Livi, and A. Suto, Phys. Rev. E 67 011102 
(2003). 

[14] R.Z. Sagdeev and G.M. Zaslavsky II Nuovo Cimento B 
97, 119 (1987). 

[15] K. Manolakou, A. Anastasiadis and L. Vlahos Astronomy 

and Astrophysics 345, 653 (1999). 
[16] Consider the case with two kinds of obstacles with re- 



spective dacay rates Ai and A2. Let Pi{t) and P2{t) 
the probability to have an encounter between times t 
and t -\- At for events of types 1 or 2 respectively. Let 
N^{t) = 1 - Pi{u)du for i = 1,2, the probabil- 
ity not to have had an encounter of type i before t. 
Then dNi/dt = ~K{Ni -\- N2). Thus d{Ni -\- N2)/dt = 
— (Ai -I- A2)(Ai -I- A2). Ni + N2 is exponential with decay 
rate Ai -I- A2 and Piit) = Aiexp(— (Ai -I- \2)t). The gener- 
alization of this result to an infinite number of rates gives 
the result 

[17] P.L. Kraprivsky and S. Redner, Phys. Rev .E 56 3822 
(1997). 

[18] At second order in v/V, the term 5\v — M{a)\r^ — 
(1 - M(a))V] can be replaced by 5(v^ - v) + 
[(1 - M(a))V].a,B(5(v« - v) + [(1 - M(a))V]4(l - 

M(a))V]j9^B bS{v^ — v) where the convention of sum- 

i j 

ming on repeated indexes is used; we compute at sec- 
ond order in V/v: A(v^,V) = A(v,V) = 2apv{l - 
V/vcos{4>) + sm^{4,)V^/2v^), and \{v) = 2apv{l + 
{V'^)/Av'^) ; always at second order in V/v, we perform 
the averages in the expression of T to obtain dtf{'v,t) = 
Dv{Af -\- -ij'v.Vvf). The term v.Vv/ stems from the fact 
that a particle more frequently collides against scatters 
with opposite velocity. Writing the equation for g{v), we 
obtain expression Q. 

[19] Using the scaling solution, we have gk ~ hk{r^,£,)/t'' , 
dtgk = A{r),(,)/t''+\ Lkvgk = ^1(77, OA'"' - Thus the ne- 
glected terms are of order l/t^ times the conserved ones, 
this is a good check of our approximations. 

[20] P. Castiglione, A. Mazzino, P. Muratore-Ginanneschi and 
A. Vulpiani Physica D 134, 75 (1999). 



4 



